/*
 * EquilibriumCondition.cpp
 *
 *  Created on: 8 Aug 2011
 *      Author: Allan
 */

#include "EquilibriumCondition.h"

// C++ includes
#include <functional>
using namespace std;
using namespace std::placeholders;

// GeoReact includes
#include "../Multiphase.h"
#include "../ReactionSystem.h"

EquilibriumCondition::EquilibriumCondition() :
multiphase(0), equilibriumReactions(0)
{}

EquilibriumCondition::EquilibriumCondition(const Multiphase& multiphase, const ReactionSystem& equilibriumReactions) :
multiphase(&multiphase), equilibriumReactions(&equilibriumReactions)
{
	// The species map object to determine the indexes of the species
	const SpeciesMap& speciesMap = multiphase.GetSpeciesMap();
	
	// The indexes of the primary (jIndexes) and equilibrium (eIndexes) species
	jIndexes = speciesMap[ equilibriumReactions.Products() ];
	eIndexes = speciesMap[ equilibriumReactions.Reactants() ];
	
	// The number of primary (Nj) and equilibrium (Ne) species
	Nj = jIndexes.size();
	Ne = eIndexes.size();
	
	// Initialise the Eigen vectors with zero entries
	alphaRow = VectorXd::Zero(Nj);
	 betaRow = VectorXd::Zero(Ne);
	   
	// Initialise the equilibrium condition function
	h = Constant(0.0);
}

void EquilibriumCondition::SetTotalComposition(const string& species, double moles)
{
	// Use the method SetComposition to setup the alphaRow vector
	SetComposition(species, moles);
	
	// The stoichiometric matrix of the equilibrium-controlled reactions
	MatrixXd vej = equilibriumReactions->StoichiometricMatrix();
		
	// Get the local primary index of the given species
	const uint j = LocalIndexPrimarySpecies(species);
	
	// If the given species is a primary species, set the betaRow vector to be the j-th column of the stoichiometric matrix vej
	if(j < Nj) betaRow = vej.col(j);
}

void EquilibriumCondition::SetComposition(const string& species, double moles)
{
	// Use the method SetCompositionAuxiliar to setup this equilibrium condition
	SetCompositionAuxiliar(species, Constant(moles));
}

void EquilibriumCondition::SetActivity(const string& species, double activity)
{
	// The global index of the given species
	const Index i = multiphase->GetSpeciesMap()[species];
	
	// Define the right-hand side function for this equilibrium condition
	auto h = [=](double T, double P, const VectorXd& n, const VectorXd& a)
		{ return n[i]/a[i] * activity; };
	
	// Use the method SetCompositionAuxiliar to setup the equilibrium condition
	SetCompositionAuxiliar(species, h);
}

void EquilibriumCondition::SetAcidity(double pH)
{
	// Use the method SetActivity to setup this equilibrium condition
	return SetActivity("H[+]", std::pow(10.0, -pH));
}

void EquilibriumCondition::SetPartialPressure(const string& gas, double Pg)
{
	// Find the gaseous phase among all the phases in the reactive system
	const Phase& gaseousPhase = *std::find_if(multiphase->GetPhases().begin(), multiphase->GetPhases().end(), 
			[](const Phase& phase) { return phase.GetName() == "Gaseous"; } );
	
	// Get a reference to the species map object of the phase system to determine the indexes of the gaseous species
	const SpeciesMap& speciesMap = multiphase->GetSpeciesMap();
	
	// The indexes of all the gaseous species
	vector<Index> igases = speciesMap[ gaseousPhase.GetSpecies() ];
	
	// Define the right-hand side function h for this equilibrium condition
	auto h = [&](double T, double P, const VectorXd& n, const VectorXd& a)
	{ 
		// The total number of moles of gaseous species
		double nT = 0.0; for(uint i = 0; i < igases.size(); ++i) nT += n[i];
		
		return nT * Pg/P;
	};
	
	// Use the method SetCompositionAuxiliar to setup the equilibrium condition
	SetCompositionAuxiliar(gas, h);
}

void EquilibriumCondition::SetChargeBalance()
{
	const vector<string> primarySpecies = equilibriumReactions->Products();
	const vector<string> equilibriumSpecies = equilibriumReactions->Reactants();
	
	for(uint j = 0; j < Nj; ++j)
		alphaRow[j] = ElectricalCharge(primarySpecies[j]);
	
	for(uint e = 0; e < Ne; ++e)
		betaRow[e] = ElectricalCharge(equilibriumSpecies[e]);
	
	h = Constant(0.0);
}

const uint EquilibriumCondition::LocalIndexPrimarySpecies(const string& species) const
{
	const vector<string> primarySpecies = equilibriumReactions->Products();
	
	const uint j = std::find(primarySpecies.begin(), primarySpecies.end(), species) - primarySpecies.begin();
	
	return j;
}
	
const uint EquilibriumCondition::LocalIndexEquilibriumSpecies(const string& species) const
{
	const vector<string> equilibriumSpecies = equilibriumReactions->Reactants();
	
	const uint e = std::find(equilibriumSpecies.begin(), equilibriumSpecies.end(), species) - equilibriumSpecies.begin();
	
	return e;
}

void EquilibriumCondition::SetCompositionAuxiliar(const string& species, const EquilibriumConditionFunction& h)
{
	// Set the equilibrium condition function (h)
	this->h = h;
	
	// Get the local primary and equilibrium indexes of the given species
	const uint j = LocalIndexPrimarySpecies(species);
	const uint e = LocalIndexEquilibriumSpecies(species);
	
	// If the given species is a primary species, set the alphaRow vector
	if(j < Nj) alphaRow[j] = 1.0;
	
	// If the given species is a equilibrium species, set the betaRow vector
	if(e < Ne)  betaRow[e] = 1.0;
}
